#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

/*
*   Nesta implementação, utilizamos a matriz de coeficientes e 4 vetores,
*   X que guarda os valores dos coeficientes na iteração k+1
*   Xant que guarda os valores dos coeficientes na iteração k
*   B que guarda os valores de B*(B dividido pelo elemento da diagonal principal)
*   row_teste que armazena os coeficientes da linha a ser testada no final
*
*   row_teste foi necessario pois alteramos a matriz incial para efetuar os calculos
*
*   optamos por nao guardar o resultado de todas as iterações para otimizar o uso
*   de memória, tendo em vista matrizes grandes.
*
*/

float *aloca_vetor(int ordem);
void preenche_row_teste(int ordem,int row,float **matriz,float *row_test);

int malloc2dfloat(float ***array, int n, int m) {
    
    /* allocate the n*m contiguous items */
    float *p = (float *)calloc(n*m,sizeof(float));
    //float *p = (float *)malloc(n*m*sizeof(float));
    if (!p) return -1;
    
    /* allocate the row pointers into the memory */
    (*array) = (float **)malloc(n*sizeof(float*));
    if (!(*array)) {
        free(p);
        return -1;
    }
    
    /* set up the pointers into the contiguous memory */
    int i;
    for (i=0; i<n; i++)
        (*array)[i] = &(p[i*m]);
    
    return 0;
}


/*
*
*  Funcao que libera memoria
*
*/
int free2dfloat(float ***array) {
    /* free the memory - the first element of the array is at the start */
    free(&((*array)[0][0]));
    
    /* free the pointers into the memory */
    free(*array);
    
    return 0;
}





/*
*   Funcao que efetua a leitura dos valores da matriz de coeficientes no arquivo
*/
int preenche_matriz(FILE *fp,int ordem,float **matriz){
    int i,j;
    float valor;
    for(i=0;i<ordem;i++){
        for(j=0;j<ordem;j++){
            fscanf(fp,"%f",&valor);
            matriz[i][j] = valor;
        }
    }
    return 0;
}

/*
*   Preenche um vetor para efetuar o teste final
*/
void preenche_row_teste(int ordem,int row,float **matriz,float *row_test){
    int i;
    for(i=0;i<ordem;i++){
        row_test[i]=matriz[row][i];
    }
    
}

/*
*   Funcao que efetua a leitura dos valores do vetor X no arquivo
*/
int preenche_B(FILE *fp,int ordem,float *B){
    int i;
    float valor;
    for(i=0;i<ordem;i++){
        fscanf(fp,"%f",&valor);
        B[i] = valor;
    }
    return 0;
}

/*
*   Faz o calculo de matriz*, X0 e B*
*/
int calc_mX(int ordem,float **matriz,float *X,float *B, float *Xant){
    int i,j;
    float div;
    //loop para dividir os membros da linha pelo valor da diagonal principal da linha
    for(i=0;i<ordem;i++){
        div = matriz[i][i];
        for(j=0;j<ordem;j++){
            matriz[i][j] /= div;
        }
        Xant[i] = B[i]/div;
    }
    return 0;
}

//Funcao que efetua o pre processamento dos dados
void prepareMatrices(float **mA, float *mB, int lines, int colums){
    int i = 0, j = 0;
    for(i=0;i<lines;i++){ 
        for(j=0;j<colums;j++){
            if(i!=j){
                mA[i][j] = -(mA[i][j]); //pela forma da equacao os elementos mudao de sinal
            }
            else{
                mA[i][j] = mB[i]; //elemento do vetor b
            }
        }
    }
}

/*
*   Funcao que calcula o valor de Xk+1 e atualiza o valor de Xk
*/
void processingMatrix(int ordem, float **main, float *X, float *Xant ){
    int i, j;
    for (i = 0; i < ordem; i++) {
        X[i] = 0;
        for (j = 0; j < ordem; j++) {
            if(i!=j){
                X[i] += main[i][j] * Xant[j];
            }else{
                X[i] += main[i][j];
            }
        }
    }
}

/*
*   Retorna o erro entre as iteracoes
*/
float relativeError(int lines, float *X, float *Xant){
    int i;
    float maxD = -3000;
    float maxN = -3000;
    float tmpd = -3000;
    float tmpn = -3000;
    for (i = 0; i < lines; i++) {
        tmpn = fabs(X[i] - Xant[i]);
        tmpd = fabs(X[i]);
        if(tmpn > maxN)
            maxN = tmpn;
        if(tmpd > maxD)
            maxD = tmpd;
    }
    
    return maxN/maxD;
}



/*
*   Main
*/

int main(int argc, char *argv[]){
    
    FILE *fp;
    int ordem, fila, iteracoes;//dados do arquivo
    int i, j;//iterador
    float tol,erro = 99,valorB;//tolerancia e erro,valorB armazena o resultado da equacao original
    float **matriz;//matriz
    float *X, *B, *Xant, *row_test;//vetores
    clock_t inicio, fim;
    double tempo;
    
    printf("\n\nCalculando...\n\n");
    
    if(!(fp = fopen(argv[1],"r"))){//abre o arquivo
        printf("Erro! Verifique se o arquivo existe!\n");
        exit(1);
    }
    
    //pega dados do arquivo
    fscanf(fp,"%d",&ordem);
    fscanf(fp,"%d",&fila);
    fscanf(fp,"%f",&tol);
    fscanf(fp,"%d",&iteracoes);
    
    //aloca matriz
    //matriz = aloca_matriz(ordem,matriz);
    //aloca_matriz(ordem,matriz);
    malloc2dfloat(&matriz, ordem, ordem);
    
    //aloca vetor para efetuar o teste final
    //row_test = aloca_vetor(ordem);
    row_test = (float *)malloc(ordem* sizeof(float));
    
    //aloca vetor de iteracoes
    //X = aloca_vetor(ordem);
    X = (float *)malloc(ordem* sizeof(float));
    
    //aloca B
    //B = aloca_vetor(ordem);
    B = (float*)malloc(ordem* sizeof(float));
    
    //aloca vetor Xant usado para calcular o erro
    Xant =(float*)malloc(ordem* sizeof(float));
    
    //preenche a matriz com os dados
    preenche_matriz(fp,ordem,matriz);
    
    //preenche o vetor de teste
    preenche_row_teste(ordem,fila,matriz,row_test);
    
    //preenche o vetor B
    preenche_B(fp,ordem,B);
    valorB = B[fila];
    
    //calcula matriz*,X0 e B*
    calc_mX(ordem,matriz,X,B, Xant);
    
    //fecha o arquivo
    fclose(fp);
    
    //loop de calculo
    inicio = clock();
    
    //efetua o pre processamento dos dados
    prepareMatrices(matriz, Xant, ordem, ordem);
    
    //iterações de calculo
    for(i=0; i <= iteracoes; i++){

          //efetua o calculo do novo valor de X
        processingMatrix(ordem, matriz, X, Xant);

        //retorna o erro calculado entre X e Xant
        erro = relativeError(ordem, X, Xant);
        
        //se o erro ainda for maior que a tolerancia, atualiza Xant
        if(erro > tol){
            for(j = 0;j < ordem;j++)
                Xant[j] = X[j];
        }else{
            break;
        }
        
    }
    fim = clock();
    
    float result=0;
    
    printf("Numero de iteracoes %d.\n",i-1);
    //calcula a aproximacao
    for(i=0;i<ordem;i++){
        result = result + X[i]*row_test[i];
    }
    //calcula o tempo utilizado nas iteracoes
    tempo = ((fim-inicio)/(CLOCKS_PER_SEC/1000));
    printf("Tempo utilzado em ms %lf.\n",tempo);
    printf("RowTest %d: => [%f] =? %f\n",fila,result,valorB);
    
    
    //libera memoria
    free(matriz);
    free(X);
    free(Xant);
    free(B);
    return 0;
    
}

